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Abstract 

The purpose of this work is simulation of magnetised plasmas in the ITER 
project framework. In this context, kinetic Vlasov-Poisson like models are 
used to simulate core turbulence in the tokamak in a toroidal geometry. 
This leads to heavy simulations because a 6D dimensional problem has to 
be solved, even if reduced to a 5D in so called gyrokinetic models. Accurate 
schemes, parallel algorithms need to be designed to bear these simulations. 
This paper describes the numerical studies to improve robustness of the 
conservative PSM scheme in the context of its development in the GYSELA 
code. In this paper, we only consider the 4D drift-kinetic model which is 
the backbone of the 5D gyrokinetic models and relevant to build a robust 
and accurate numerical method. 
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1. Introduction 

The ITER device is a tokamak designed to study controlled thermonu- 
clear fusion. Roughly speaking, it is a toroidal vessel containing a magne- 
tised plasma where fusion reactions occur. The plasma is kept out of the ves- 
sel walls by a magnetic field which lines have a specific helicoidal geometry. 
However, turbulence develops in the plasma and leads to thermal transport 
which decreases the confinement efficiency and thus needs a careful study. 
Plasma is constituted of ions and electrons, which motion is induced by the 
magnetic field. The characteristic mean free path is high, even compared 
with the vessel size, therefore a kinetic description of particles is required, see 
Dimits Then the full 6D Vlasov-Poisson model should be used for both 
ions and electrons to properly describe the plasma evolution. However, the 
plasma flow in presence of a strong magnetic field has characteristics that 
allow some physical assumptions to reduce the model. First, the Larmor ra- 
dius, i.e. the radius of the cyclotronic motion of particles around magnetic 
field lines, can be considered as small compared with the tokamak size and 
the gyration frequency very fast compared to the plasma frequency. Thus 
this motion can be averaged (gyro-average) becoming the so-called guid- 
ing center motion. As a consequence, 6D Vlasov-Poisson model is reduced 
to a 5D gyrokinetic model by averaging equations in such a way the 6D 
toroidal coordinate system {r,0, (j),v^\,v±, a) becomes a 5D coordinate sys- 
tem (r, 0, (p, , /i), with I'll the parallel and v± the perpendicular to the field 
lines components of the particles velocity, a the angular velocity around the 
field lines and n = m v\/2B the magnetic momentum depending on the 
velocity norm |f_L|, on the magnetic field magnitude B and on the particles 
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mass m which is an adiabatic invariant. Moreover, the magnetic field is 
assumed to be steady and the mass of electrons me is very small compared 
to the mass of ions mj. Thus the cyclotron frequency Wj^g = <li,e B/nii^f. is 
much faster for electrons than for ions >> Wj. Therefore the electrons 
are assumed to be at equilibrium, i.e. the effect of the electrons cyclotronic 
motion is neglected and their distribution is then supposed to be constant 
in time. The 5D gyrokinetic model then reduces to a Vlasov like equation 
for ions guiding center motion: 



8t +^-U-f"j+««.U^'-j='' w 

where /^(X, Wji) is the ion distribution function for a given adiabatic in- 
variant /V, with X = (r, 9, 0), velocities dX/dt and dv^^/dt define the guiding 
center trajectories. UV(^x,v^^)'{d^/dt,dv\^/dty = 0, then the model is termed 
as conservative and is equivalent to a Vlasov equation in its advective form: 

t + f ■^-^' + >«.(^-)-<'- 

This equation for ions is coupled with a quasi-neutrality equation for the 

electric potential on real particles position, with R = X — (with 
the Larmor radius) : 

1 e 
-V± • (noV±$) + —($-<$ >0,4,) = J f^,dpd.v\\ - no (3) 



where no is an equilibrium electronic density, Tg the electronic temperature, 
e the electronic charge, k the Boltzmann constant for electrons and coi the 
cyclotronic frequency for ions. 

These equations are of a simple form, but they have to be solved very ef- 
ficiently because of the 5D space and the large characteristic time scales 
considered. This work is then a contribution in this direction, following 



Grandgirard et al who develops the GYSELA 5D code that solves this 5-D 
gyrokinetic model, see [5] and [6]. Looking at the model, one notices that 
the adiabatic invariant /i acts as a parameter. Therefore for each ^ we have 
to solve a 4D advection equation as accurately as possible but also taking 
special care on mass and energy conservation, especially in this context of 
large characteristic time scales. The maximum principle that exists at the 
continuous level for the Vlasov equation should also be carefully studied at 
discrete level: 

min(/(x„r)) < < max(/(xi,r)) 

with f(xi,t^) the value at Xi in cell i at time t". There is no physi- 
cal dissipation process in the gyrokinetic model ([T]) that might dissipate 
over /undershoots created by the scheme and the loss of this bounding ex- 
trema of the solution at t"'^^ > may even eventually crash a simulation. 
Those studies will be achieved in this paper on a relevant reduced model, 
the 4D drift-kinetic model described in section |4j which has the same struc- 
ture than equations ([T]) . The geometrical assumptions of this model for ion 
plasma turbulence are a cylindrical geometry with coordinates {r,6,z,v^^) 
and a constant magnetic field B = Bz e^, where ez is the unit vector in 
z direction. This 4D model is conservative and will be discretized using a 
conservative semi-Lagrangian scheme, the Parabolic Spline Method scheme 
(PSM, see Zerroukat et al [12j and [13j ). It is a fourth order scheme which 
is equivalent for linear advections to the Backward Semi-Lagrangian scheme 
(BSL) currently used in the GYSELA code (see Grandgirard et al \&\) and 
introduced by Gheng-Knorr [2] and Sonnendriicker et al [TO]). This conser- 
vative PSM scheme based on the conservative form of the Vlasov equation 
will be described in section [4] and properly allows a directional splitting. 

4 



In this paper, the BSL and PSM schemes wih be detailed with an empha- 
sis on their similarities and differences. We will see that one difference is 
about the maximmn principle. The BSL scheme satisfies it only with a 
condition on the distribution function reconstruction and the conservative 
PSM scheme does not satisfy it without an extra condition on the volumes 
conservation in the phase space. The last condition is equivalent to try 
to impose that the velocity field is divergence free at the discrete level. A 
scheme is given to satisfy this constraint in the form of an equivalent Finite 
Volume scheme. Moreover, we have designed a slope limiting procedure, 
Slope Limited Splines (SLS), to get closer to a maximum principle for the 
discrete solution, by at least diminish the spurious oscillations appearing 
when strong gradients exist in the distribution function profile. 
The outline of this paper is the following : in section [2] will be recalled some 
important properties of Vlasov equations at the continuous level. Then BSL 
and PSM schemes will be described and compared, according to properties 
of the discrete solutions. In section |3| a numerical method will be given to 
improve the respect of the maximum principle by Vlasov discrete solutions 
when using the PSM scheme and particularly to keep constant the volume 
in the phase space. In section |4j practical aspects of the PSM scheme use 
will be described in the context of the 4D drift-kinetic model and at last we 
will comment on numerical results. 
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2. Semi-Lagrangian schemes for Vlasov equation 

2. 1 . Basics of the Vlasov equation 

Let us consider an advection equation of a positive scalar function /(x, t) 
with an arbitrary divergence free velocity field: 

at/ + a-V.(/) = 0, 
V • a = and f{x,t) > 

with position x E and a{x, t) G R-^ the advection velocity field. 
The solutions satisfy the maximum principle: 

0</(x,t) <max(/(x,to)) (5) 

X 

for any initial time to < t- 

Since V • a = 0, we can also use an equivalent conservative formulation of 
the Vlasov equation: 

dtf + V,. • (a /) = 0, (6) 

For more details, see Sonnendriicker Lecture Notes [TT]. One obvious prop- 
erty of this conservation law (Reynolds transport theorem) is to conserve 
the mass in a Lagrangian volume Vol{t), by integrating the distribution 
function on each Lagrangian volume element dO: 

dtm = dt I f{x,t)dQ = 0. 

Jvol{t) 

Let us introduce the convective derivative dt{.) = dt{.) + a ■ V x{-)-, thus ([g]) 
becomes: 

dtf + / V, • a = 0. (7) 

Considering a Lagrangian motion of an infinitely small volume Vol{t), we 
have dtm = dt{f Vol) = 0, thus we obtain: 

dtVol 



Vol 
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• a. (8) 



Obviously, a divergence free flow V^^ -0 = conserves a Lagrangian volume 
in its motion. 

2.2. Maximum principle in the BSL and PSM schemes 
2.2.1. Backward semi- Lagrangian (BSL) 

Let us consider a Vlasov equation in its non conservative form: 

dtf + a- V,./ = 0, (9) 

with /(x, t) a scalar function, position x G and a{x, t) G the ad- 
vection field. The BSL scheme, see Sonnendriicker et al [.lOj, is based on 
the invariance property of function / along characteristic curves to obtain 
values Z""*"^ at time t""^^ from the values /" at t": 

jn+l t"+l)) = /" (X(x"+\ t")) , (10) 

with x the Eulerian coordinates and the characteristic curves X defined as 

dX(x,t) , . 

-^-^=a{x,t) (11) 

with the initial position x = X{x,t^) at t". Let us locate the discrete func- 
tion values /f = r (^«^\i"+^)) at mesh nodes = X(x^+\t"+i). 
We solve the following nonlinear system which is a second order approxima- 
tion of dtX{t) = a{x,t): 

= {Xix^+\t'^+^)+X{x^+\e^)) /2, 
X{x^+\n = X(x^+\ t"+i) - At a (Xf , (12) 

r+i (X«+\r+i)) = /« {Xix^+\n) , 

with At = t""*"^ — t" . The function fj^{x) is a reconstruction of the solution 
/"(x) according known values at nodes x"^^ using cubic splines basis func- 
tions on the domain to obtain the value at x" = X{x^^^,t"'), which is not 
a mesh node in general. 
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Properties of the BSL scheme. This scheme is formany fourth order in space. 
It is second order in time using for instance a Leap- Frog, Predictor-Corrector 
or Runge-Kutta time integration. Mass is not conserved by this scheme, 
because it has no conservative form. However, an approximated maximum 
principle is satisfied. Let us consider fj^{x) the cubic sphne interpolation of 
the distribution function f{x,t'^) at time t", we have for any x: 

r+i (x(x,r+i)) = 4"(x(x,t")). 

It then naturally appears a "discrete" maximum principle: 

min(4"(x,r)) < r+\x) < max(/,"(x,P)). (13) 

X X 

Comparing with the property (pi), we have here minf/I'fx, t")) 7^ and 

' — ' X 

max(/^(x, t")) 7^ max(/"(j;, t"")) because the cubic spline reconstruction 

X X 

does not satisfy a maximum principle. If we have a manner to enforce this 
property to this reconstruction, a maximum principle is granted for the BSL 
scheme. No directional splitting is allowed since the BSL scheme is based 
on the non conservative form of the Vlasov equation, see 

2.2.2. Semi-Lagrangian Parabolic Spline Method (PSM) 

Let us consider a Vlasov equation in its conservative form: 

dtf + Vx-{af)= 0, (14) 

with f{x,t) a scalar function, position x € and a{x,t) G MP the advec- 
tion field. Notice that with the hypothesis • (a) = 0, conservative form 
(14) and non-conservative form ^ of the Vlasov equation are equivalent. 



The PSM scheme, see Zerroukat et al [12] and |13j, is based on the mass 
conservation property of function / in a Lagrangian volume Vol to obtain 
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the value /"+^ at time 

[ f{x,r+^)dn= [ f{x,e)dn, (15) 

Jvol"+^ Jvol" 

dX(x, t) 

with the characteristic curves X defined as — ^ — = a(x, t) and = 

X{x,t'^) with the initial time, and the volume Fo/" = t") such 

that t"^"*^) G yoZ""^"*^} defined by the Lagrangian motion with the 

field a{x,t). The important point is that this conservative formalism prop- 
erly allows a directional splitting without loosing the mass conservation. 



Indeed, equation (14) may be solved with D successive ID advections still 



of conservative form: 

dtf + d^,{akf) = 0,ke[l,D]. (16) 

We then approximate a ID equation for each direction k using the conser- 
vation property. Omitting subscript k, the PSM scheme writes in ID as 
follows: 



fix,e+')dx= I f{x,ndx, (17) 

n+1 ' 



with x^^y^ = X(x"^jy2' *"^^) settled as the ID mesh nodes and x^_^_^^2 = 
^(^7+1/2'' associated foot of the characteristic curve, Vol^ = [r" 



n+1 1 
i-1/2' •^i+l/2i- 



and Vol"+^ = [x"+\o,2; 



Let us define the unknowns of the scheme as the average of / in cell i 



-n+1 1 r«+i/2 



i-l/2 



and the primitive function 

F"(z)= / f{y,ndy, (19) 



1/2 

with the uniform space step Ax = x^^^^^ — x^^^^^ and Xi/2 an arbitrary 
reference point of the domain and for instance the first node of the grid 
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{xj_i/2}j=i,Af+i- Therefore, one has to solve a nonhnear system, which is 



similar to the BSL one, to obtain a discrete solution of equation (17) that 
writes: 

X{x^},„n = - At a (x-:'/lt-+yA , (20) 



with the time step At = t'^'^^ — t" and the uniform space step Ax = x^^y^ — 

^i-l/2 ■ 

The computation of the reconstructed primitive function FJ^{x) is based on 
values at mesh nodes x^^^^^' 

i 

k=l 

Then this set of values is interpolated by cubic splines functions to obtain 
an approximated value FJ^{z) of the primitive function F^{z) at any point 
z of the domain: 



FJi{z)^F^{z)= / f{y,ndy. (21) 

•^3:^1/2 

Properties of the PSM scheme. This scheme is formally fourth order in space 
and strictly equivalent to the BSL scheme for constant linear advection, see 
[3]. It is second order in time using for instance a Leap- Frog, Predictor- 
Corrector or Runge-Kutta time integration scheme. Mass is exactly con- 
served by this scheme for each ID step k of the directional splitting. How- 
ever, no maximum principle does exist for each step k even for the exact 
solution: in general dx^a^ 7^ 0, even if the velocity field is divergence free 
V • a = 0. Let us consider the scheme in D dimensions of space: 

dtf + V- (a /) = 0, (22) 
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with f{x,t) a scalar function, position x G and a{x,t) G the advec- 
tion field. Let us consider a cell i, where the solution is described at time 

by its average in cell i, the PSM scheme then writes: 

f^+Woi^+^=[ fix,e+')dn= [ f{x,ndn. (23) 

We thus obtain the following relation: 

yn+l ^ yn* Volj , . 

with the average of the distribution function in the Lagrangian volume at 
time t": 

7: 



' Vol, 



Oh J Vol? 



Here clearly appears two conditions, both difficult to satisfy especially in 
the context of a directional splitting, to have a maximum principle defined 
as follows: 

min(7;) < < max(7;). (25) 

1. Maximum principle on the distribution function in Vol": 

min(7;)<7r<max(7;). (26) 

2. Conservation of volumes in the phase space at the discrete level: 

Volf = Vol^+^. (27) 

The first condition is difficult to ensure in general, because a maximum prin- 
ciple should be satisfied for any average of the distribution function on an 
arbitrary volume Volf. Moreover, in the context of a directional splitting, 
it is impossible to satisfy a maximum principle for a ID step k, because it 
does not exist at the continuous level since in general dkak 7^ 0. Therefore it 
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is probably impossible to recover a maximum principle of the reconstruction 
after all steps of the directional splitting. 

The second condition is true at the continuous level while V • a = 0, since 



we have dtVol = Vol V • a = 0, see equation ([8j) in section 2.1 As well as 
for the first condition, in the context of a directional splitting, it is difficult 
to ensure a constant volume evolution Vol'^ = Vol^^^ after all steps of the 
directional splitting, where compressions or expansions of the Lagrangian 
volume occur successively. 

As a consequence, we will propose a form of the conservative PSM scheme 
that does not use a directional splitting. However we will not write the Semi- 
Lagrangian form of the PSM scheme in D dimensions of space because it 
is costly in computational time, because of the reconstruction step, and 
it is difficult to handle with arbitrary coordinate systems. The solution 
we choose is to use an equivalent Finite Volume form of the PSM scheme 



described in section 3.2, which is locally ID at each face of the mesh. It is 



therefore possible to design ID numerical limiters to try to better satisfy the 



maximum principle condition (26). Moreover, we will show that this form 



allows an exact conservation of the volumes in the phase space (27). The 
maximum principle and therefore the robustness of this scheme will thus be 
considerably improved. 

3. Maximum principle for the PSM scheme 

3.1. Numerical limiters for the distribution function reconstruction 

Enforcing the first condition on the maximum principle of the distribu- 



tion function reconstruction (26) can be really costly in computational time. 
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Instead of trying to correct the cubic spline reconstruction, we will reduce 
the spurious oscillations, generated by high order schemes when strong gra- 
dients appear in the distribution function profile, by using a classical Van 
Leer like slope limiting procedure, see for instance LeVeque [9j. We propose 
here to measure the gradients in the flow and to add diffusion where they 
are detected. The diffusion is added by mixing the high order PSM flux 
with a first order upwind flux. The evaluation of the gradient is given by 
the classical function 9 and we estimate the diffusion needed with a function 
7(0) G [0,1] based on a minmod like limiter function (see Fig. [T]). The 
resulting limiter we propose here is called SLS (Slope Limited Splines), see 
[7] for details: 

= 7(^.+l/2) Cl/2 + (1 - 7(^.+l/2)) <^ZT/f 



where 



iUpwind 
J+1/2 



/Tl pTL pTl nil 

i + J i+1 ■ ( \ J i+1 ~ Ji 



«i+i/2 7) sign(ai+i/2 



We define as the classical slope ratio of the distribution which depends 

on the direction of the displacement: 



^i+1/2 



if , 1 /9 > 



/I jil L 

i / i—1 

=n =n II O.i+1/2 

J i+1 ~ Ji 

Ji+2 ~ J i+1 -f , n. 

^=n =r!r II «i+l/2 < ^ 

J i+1 ~ Ji 



However, the classical limiter minmod, where 7i+i/2 = max(0, mm{9i_^^i/2, 1)), 
set 7 to when ^ < 0. That means that the scheme turns to order 1 when 
an extrema exists, i.e. the slope ratio 6 < 0. These extrema are thus quickly 
diffused and that leads to loose the benefits of a high order method. For 
SLS, the choice is to let the high-order scheme deal with the extrema and 
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Figure 1: 7 function for the SLS limiter 

only add diffusion when strong gradients occurs, i.e. the slope ratio ^ ~ 0. 
We also introduce a constant K in relation to control the maximum slope 
allowed without adding diffusion, i.e. mixing with the upwind scheme, see 
figure [Tj 

7i+i/2 = max(0,min(K|6'i+i/2|, 1)). (28) 
with the constant K = 5 experimentally settled. 

We present in Fig. [2] the results of the linear advection of a step function 
with the standard PSM scheme with and without the SLS limiter (K=5). 
The domain is meshed with 70 cells with periodic boundary conditions and 
the displacement is set to 0.2 cell per iteration. One can see that as any high 
order scheme, the PSM scheme produces spurious oscillations at the discon- 
tinuity or at a stiff gradient location. The SLS limiter do well with K = 5 
to reduce these oscillations without introducing much diffusion. However, 
a maximum principle is not granted. This limiter has been further studied 
and compared with other limiters in the report [7] . 
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Figure 2: Linear advection of a step function. The exact solution is green and the 
numerical results are in blue, above the standard PSM scheme and below the PSM scheme 
with SLS limiter with AT = 5. 

3.2. Finite Volume form of the PSM scheme 

Let us consider the ID conservative advection equation of the form: 

dtf + d^{a^ f) = 0, (29) 

with f{x, t) a scalar function, position x G M and a^^x, t) G M the advection 
field. We recall that X{x^^^,^,t"'^^) = Xi^i/2 is the position of the mesh 



"foot" 



node i + 1/2. Let us set the notation X{x^^^^^,t'^ 
position at on the characteristic curve. Let us rewrite the PSM scheme 



(20): 



/rAx = F,"(x*_,i/2)-KK-i/2)- 
with the primitive function FJ^{z) is interpolated by cubic splines at mesh 
nodes such that FJ^{z) ^ F^{z) = J^^^^ f^{y)dy at fourth order in space. 
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Let us make appear explicitly in ID the fluxes at cell faces i it 1/2, by 
introducing the primitive values at cell faces F^{xi±i/2)- 

/r'Ax = [fI:{x*^,/,) - - (Fr«-i/2) - K(^^-i/2)) +7>x 

with /"Arc = F^"(x,+i/2) - FJ^{x,_y2)- It yields 

/r' - 1: , (^.+1/2) - i^/r«+i/2) (^.-1/2) - (^*-i/2) _ Q 

At AxAt AxAt 

(30) 

The PSM fluxes at cell faces z it 1/2 clearly appear: 

-fn+l -fn (rPSM _ ^PSM 

l^^A + ^l±lIl^lzlIl=o (31) 
At Ax ^ ^ 

with 

•^^'1+1/2 

A simple Taylor expansion shows that this PSM flux, which consists in a 
cubic spline approximation of the integral of /"(x) along the characteristic 
curves at cell faces, is a consistent approximation at node i + 1/2 of the 
continuous flux <1> = Uxf in equation (29), i.e. ^f^if2 ~ {(^xf)i+i/2- 
Moreover, this flux is an approximation of the integral of {a^f^ on cell faces. 



Coming back to (29) and integrating on the cell volume Vol^~^^ = A^Ax, 



with Ax the bounding faces Tj_-^i/2 area transversal to x of cell i: 

+ [ dxiax f)dn = 0. (33) 

JVol"'+^ 



We obtain using Green formula: 



df . 1 



+ / /(ax • nx)dT = 0. (34) 
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with dVol^~^^ = rj_i/2 U ^4+1/2 the surface bounding Vol"^^ and Hx its 
outgoing normal. 



dt Voir' \Jr' Jr' / 



(35) 



Comparing with formula (31 ), we see that the PSM flux is an approximation 



of the flux average at cell faces: 



because Vol'l^^ = A^Ax 



The extension to D dimensions of space is then straightforward, because 
it only consists in adding the fluxes through the faces of a cell in every 
direction d. Considering the Vlasov equation in its conservative form in 
dimension D and for an arbitrary coordinate system: 

^+V.(Ja/)=0 

with J the geometric Jacobian of the cell, f{x, t) a scalar function, position 
X G and a{x, t) G the advection field. In a classical way in Finite 
Volume methods, we integrate this local equation on the cell Cj of volume 
Vol^^ and we use the Green formula: 

/ %Jdx + I {a- nd)fJdx = 0, (37) 
JVoP+i ot Jvaec^ 

with Frf the face of cell i perpendicular to direction d of area and of out- 



going normal unit vector n^. Using the ID flux formula (36) in direction 
and a first order time discretisation, we obtain the following Finite Volume 
scheme: 

-jn+l -j-n 
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with 



fl = vr^l ,f{y,nJ{y)dy. (39) 



Vol' 



n+l 



Vol]' 



and (^PS^ the flux that goes through Vd'- 

= AdJr ^'^'^^y^ * )J{y)dy, (40) 

with x*^ = Xd — At {a{xd) ■ rid) the foot of the ID characteristic curve and 
F^^(z) the primitive function at t" reconstructed using cubic sphnes in the 
direction of Ud'- 

FhA')^F2i^)= f{y,nJiy)dy. (41) 

''Xd,l/2 

We see here a Finite Volume form of the Semi-Lagrangian PSM conservative 
scheme. This equivalence is however restricted, because this Finite Volume 
form is submitted to a CFL condition as any scheme of this form: 



At < min 



Axd 



d \ max(a2(xrf)) 



Moreover, as we will see in section 3.3, the Lagrangian volume evolution is 



here approximated by the cell faces motion only in their normal direction, 
instead of the general motion as it is described in the Semi-Lagrangian 



formalism (15). It is the same volume evolution as the Semi-Lagrangian 
method with ID directional splitting. However it is the classical Finite 
Volume formalism and it is the key point that will permit to enforce a 
divergence free evolution of the flow. 



Notice that the Finite Volume form (38) can be directionally split keeping 



exactly the same result. Indeed, it only consists in adding the flux in two 
successive operations instead of in one. As an example, let us consider the 
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2D (x, y) case of a cartesian mesh: 



r J- n" 

''^'i +A+l/2,j*i+l/2,i + A-l/2j'Pi-l/2,i (42) 



I 4^ ^P5A/ , AV (hPSM 

^^i,i+l/2*iJ+l/2 ^ ^i,j-l/2*i,i-l/2 



0. 



It is stricly equivalent to use the directional directional splitting: 

T/„/"+l f i ~ fi _i_AX S)^SM I /IX fl)PSM _ p, 



At ' -^i+l/2,j^i+l/2,i A-l/2,i^j-l/2,i - 

+1 -jnx (,4dj 

'/O'i +^i,i + l/2*ij + l/2 + ^iJ-1/2%-1/2 - U- 

with the only condition that all fluxes ^^^'^^ are computed using /" and a^' 



at time as in the unsplit scheme ( 42 ) . 



3.3. Conservation of volumes in the phase space 

The second condition to have a maximum principle for the PSM scheme 



is to satisfy the multi-dimensional condition (27) of conservation of volumes, 
i.e. Vol" = Vol""^^. Equation Q showed that at continuous level the 
volume is constant in its evolution in the phase space if the advection field 
is divergence free. Therefore we will study the PSM scheme to find out 
a divergence free condition that should be satisfied at the discrete level 
V'^ • a = in such a way Vol'"' = Vol"~^^, in the same way V • a = at 
the continuous level. With the idea of making appear the total evolution 
of volumes between and t""^^, we will use the Finite Volume form of the 



PSM scheme given in section 3.2 for the 2D polar coordinate system. Let us 
consider radial r and orthoradial 9 directions, with constant space steps Ar, 
A6 and the volume of the cells Volij = r-iArAO with the mean radius 
of the cell. We consider that the mesh in polar coordinates has locally 
no curvature, i.e. each mesh is a trapezium with straight edges. This is 
important to be noticed to write the Finite Volume scheme and calculate 
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Figure 3: Polar cells in black and the Finite Volume representation of the cells in red as a 
trapezium with the volumes swept by the cell faces in their normal motion. The direction 
for all faces motion is drawn outward, but it could be inward as well, function of the 
velocity field. 

the volume swept by the cell edges, see Fig. [3j Let us set a velocity field 
{ar{r, 9), r ag{r, 9)) such that: 

Vr,6» • a = -dr{r ar) H — dg{r ag) = 0. 

Let us write the conservative advection equation in polar coordinates: 

dt{rf) + dr{r ar f) + dg{r ag f) = 0. (44) 

Notice that the geometric Jacobian J = r for polar coordinates. 

The PSM scheme without directional splitting in the Finite Volume form 



(38) reads here: 



Vol 



-jn+l -jn 

J iJ ~ J iJ ,Ar (f,PSM,r _ .r ^PSM,r 

»J +^i+l/2,j^»+l/2,i ^i-l/2,i^i-l/2,j (45) 

JJ+1/2 iJ+1/2 «,J+l/2 J,j-l/2 
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with ^^'^ *^afi/2'j positioned at cell faces center and with cell i,j of 

volume Volij = ViArAO and faces areas A^^^^i^^ = AO and ^^^-1-1/2 ~ 
The cell averaged values of / used in the scheme are: 

J =-^[ f{r,e,t)rdrde. (46) 
Using the integral form (40) of the fluxes: 



r,ArA0(/"+' - f\ 



+Ae / f{r,ei,e)rdr - AO / f{r,6i,r)rdr 

+Ar f{n,e,nnde- Ar f{n,e,nnde = Q. 



(47) 



J + 1/2 J-1/2 



Let us introduce the volumes swept by each cell face in its normal motion 
in accordance with the Green formula and the way of computation of feet of 
characteristic curves normal to cell faces, i.e. without taking into account 
the tangential motion at the cell faces or the curvature of the mesh, see Fig. 

13 

SVoll^,/2 = n At (^^.±1/2 - e;±i/2). 
Therefore we obtain: 

/ /f +^ rdrde = [ f{r, 6, T) rdrdO 

with 

We here recover a discrete mass conservation formulation. To obtain Vokj 
Vol'^j, and thus preserve a constant function, it yields: 
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(48) 



Using 6Vol^ definitions, we thus obtain a discrete divergence formulation in 
polar coordinates • a = to be nullified: 

1 ri+i/2ar{ri+i/2,Sj) - ri^i/2ariri-i/2,0j) ^ 

r ^ r ^ (49) 

1 riae{r„ 6*^+1/2) - r,ae{r„ dj-1/2) _ ^ ^ 

n " ■ 

with the following first order definition of the characteristic curves feet com- 
putation: 

ri±i/2 - A±\I2 = «r(r-i±i/2, ^j) and 6'j±i/2 - 0*^i/2 = (^eiu, 9j±i/2)- 

As a conclusion, we have presented a general methodology for any co- 
ordinate system to compute the associated discrete divergence free condi- 
tion, by using the approximation of cell edges by straight lines and by only 
considering the normal to cell faces motion of the volume as it has to be 
when invoking the Green formula in this Finite Volume framework. The 



discrete divergence formulation (49) is independent of the time integration 
method. It is a discrete consistent relation for the advection field of the 
form V'^ -0 = 0. It should be satisfied to get the conservation condition on 
volumes Vol^ = Vol^~^^ in the phase space, which is necessary to obtain a 
maximum principle for the PSM scheme or actually for any Finite Volume 
scheme. This condition is also necessary when using the Semi-Lagrangian 



PSM scheme with directional splitting as described in section 2.2.2 as well 



as when using the Finite Volume form described in section |3.2[ In Fig. H] 



we compare the results of a 4D drift-kinetic benchmark (see section 4.3.1 



for details) obtained with the Semi-Lagrangian PSM scheme 2.2.2 with an 



advection field computed: first in such a way the discrete divergence free 



condition (49) is satisfied and second with an advection field computed by 



cubic spline interpolation without satisfying this condition. 
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Figure 4: Result at time t = 60 with (left) the advection field computed in a way 
(see section |4.3.1[ ) that satisfy the discrete divergence condition | |49[ ) and with (right) 



the advection field computed with cubic splines, which do not satisfy this condition ( 49 1 



Respecting condition ( 49 1 for the advection field not only leads to a better respect of the 



maximum principle, it is actually necessary to ensure the stability of the scheme. The 
result in figure |4] diverges from realistic physics. 



4. Use of the PSM scheme in a 4D drift-kinetic code 

4-1. Drift-kinetic model 

This work follows those of Grandgirard et al in the GYSELA code, see 
[5] and [6] . The geometrical assumptions of this model for ion plasma tur- 
bulence are a cylindrical geometry with 4D coordinates {r,9,z,v^\) and a 
constant magnetic field B = Bz ez, where is the unit vector in z direc- 
tion. The model is the 4D Drift-Kinetic equations described inGrandgirard 
et al [6j: 

dr dO dz dvn 

di = ''Tt = = "II' ^ = m,''' 

with vgc = {E X B)/B'^ and E = — V<I> with <I> the electric potential. 
The 4D Vlasov equation governing this system, where the ion distribution 
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function is /(r, ^, z, ,t), is the following: 

dtf + VGcA-f + VGCedef + v\\dj + ^E,d, f = 0. (51) 

This equation is coupled with a quasi-neutrality equation for the electric 
potential ^{r,6,z) that reads: 

with rii = f{r,9,z,v\\)dv\\ and constant in time physical parameters no, 
Qq, Te and e. Let us notice that the 4D velocity field a = (fcc^i ^^GCg) ^||) 
q/rrii E^Y is divergence free: 

V • a = ^9r(r- VGcJ + ^deivGCg) + dzV\i + dy^^ {q/nii E^) = (53) 

because of variable independence dy^^Ez = c?„|| {dz^{r, 9, z)) = and dzV^^ = 
and we have vgc = {E x B)/B'^, with E = — V<I> and B = Bz e^, thus 



VGCr = and vcce = ^ (^r^) (54) 



and ^ ^ 

V, - a = -5r(r VGCr) + -de{vGCg) 

^ 1 ^ (55) 

= ^ (-l/r)ae$) + 9e (9,«>)) = 0. 

Therefore, one can write an equivalent conservative equation to the preced- 



ing Vlasov equation (51): 



dtf + dr{vGCr f) + deivGCo f) + 5.(^-11 /) + —Ez / = (56) 

4-2. Computation of a divergence free velocity field at the discrete level 
We have obtained a discrete form of the velocity field divergence to 



nullify (49), as a necessary condition to obtain a numerical solution with a 
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maximum principle. We saw in (53) that V • a = is satisfied equivalently 
if VrO ■ a = d ( [SS] ) is satisfied and this is stih true at the discrete level 
(independence of variables). Therefore, the velocity field should nullify the 



discrete polar divergence ( 49 ) : 



1 ri+i/2ar{ri+i/2,Gj) - ri_ii2ar{ri.ii2,0j) ^ 

ri Ar 

1 riae(ri,6'j+i/2) - riaB{ri,9j_ii2) _ 



(57) 



n AO 

with 

Or = dr/dt = VGCr = — Irde^ and ae = d6/dt = vcce/r = —^dr^, 



using definitions given in (54) 



Proposition 1. Let us define the electric potential at the nodes of the mesh 
^j+i/2j+i/2; whatever the way it is computed. Let us set the following nat- 
ural finite difference approximation for the velocity field: 

, a ^ _ -1 ^»+l/2,j+l/2 - ^i+l/2,j-l/2 

0'r[rij^ll2,0j} - - - — 

^1+1/2 AO ^gg^ 

, . 1 ^i+l/2J+l/2 - ^»-l/2J+l/2 



With this approximated velocity field, the approximation of^^e -0 = (51) 
is satisfied. 



The proof is easy, we just have to put the velocity field (58) in (57) to see 
that all terms annulate each others. 

Remark 2. Notice that the electric potential $ should he computed at nodes 
{i lb 1/2, j lb 1/2) of the mesh to obtain velocities at the center of cell faces 
{i lb 1/2, j) and {i,j ±1/2). It is well adapted to the PSM schemes, where 
the displacement should he calculated at cell faces. 
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4-3. Numerical tests 

4-3.1. Drift-kinetic 4D model, PSM schemes comparison 

In this section, we will compare the numerical methods on a 4D drift- 
kinetic benchmark, following the paper Grangirard et al [6j. The model is 



described in section 4.1 We will compute the growth of a 4D unstable tur- 
bulent mode. The benchmark consists of exciting the plasma mode (m,n), 
with m the poloidal mode (O) and n the toroidal mode (z). The initial 
distribution function is the sum of an equilibrium and a perturbation dis- 
tribution function / = /eg + 5f. The equilibrium distribution function has 
the following form: 

and the perturbation 6f 



6f{r,0,z,v\i) = fegir,v\i) g{r) h{v\\) 6p{6,z) (60) 

with g(r) and /i(f||) two exponential functions and 

dp{9, z) = e cos ( — — z + mO 

with Lz the length of the domain in z direction, rrii, Ti(r), no(r) physical 
constant profiles, see [6j for details. We have set here m = 16 and n = 8. 

4.3.2. Algorithm 

At the beginning of the time step, the distribution function /(x,t;||,t") 
is known at time t", with x = {r,6,z). The time step At = t""*"^ — t" is 
computed at each step with the CFL like condition: 

M = min I CFLd ^^'^ 
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with the coefficients CFL^ = CFLq = 0.5, because the flow is highly non- 
linear in (r, 0) planes thus characteristics should not cross each others during 
one time step, and CFLz = CFL^^^ = 8 because it is linear advection in 
direction z and f y so characteristics can not cross each others and then we 
allow a maximum displacement of 8 cells. Actually excluding the linear 
phase, the most restrictive directions for the time step are r and 9, in such 
a way this last value (8) has a minor importance compare to the leading 
parameters CFLy = CFLq = 0.5. 

The operator splitting between the quasi-neutral equation and the Vlasov 
transport equation is made second order using a Predictor-Corrector scheme 
in time: 

1. Time step At computation. 

2. Quasi- neutral equation ( |52[ ) solving at using the distribution func- 
tion /" (actually the density) to obtain the electric potential <I>"(x) 
at time t"^'-'^/^. The advection field a(x,t") is computed with <I>"(x) 



according to equation ( 53 ) and using formula ( 58 ) . 

3. 4D Vlasov equation solving at with time step At/2 to obtain the 
distribution function /"^"^/^(x) at time using the advection field 
a(x,t"). 

4. Quasi-neutral equation (52) solving at using the distribution 
function (actually the density) to obtain the electric potential 
^n+i/2^^-j |.jj^^g The advection field a{x,t'^^^/'^) is computed 
with <I>"+-'^/^(x) according to equation (53) and using formula (58) . 



5. 4D Vlasov equation solving at with time step At to obtain the 
distribution function /"+^(x) at time t""*"^ using the advection field 
a(x,t"+V2). 
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In the two following paragraphs, we describe the schemes for the 4D 



Vlasov equation (56) solving of the algorithm with At* = At/2 in the pre- 



diction step and At* = At in the correction step. 

4D Semi-Lagrangian PSM sheme with directional splitting. 

• PSM ID advection of /(x,U||,t"') in direction with velocity 0"^ and 
time step At*/2 to obtain /(x, , t""/^). 

• PSM ID advection of f{x, f y , t'^H''^) in direction z with velocity a" and 
time step At*/2 to obtain /(x, , t^/^). 

• PSM ID advection of f{x,v^i^,t^/'^) in direction 9 with velocity and 
time step At*/2 to obtain /(x, , t^/^). 

• PSM ID advection of /(x, , t^/'^) in direction r with velocity a" and 
time step At* to obtain /(x, , f"). 

• PSM ID advection of /(xjfujf) in direction 6 with velocity and 
time step At*/2 to obtain /(x,W||,t^). 

• PSM ID advection of /(x,V||,t^) in direction z with velocity a" and 
time step At*/2 to obtain /(x,V|j,t^). 

• PSM ID advection of /(x,f||,t^) in direction vy with velocity a^^^ and 
time step At*/2 to obtain /(x, wy , ) = /(x, wy , t"+i). 

Each PSM ID advection is achieved using the standard ID semi-Lagrangian 



PSM scheme as described in section 2.2.2 The directional splitting is second 



order by using a Strang like decomposition. Since we use here a directional 
splitting and a second order scheme in time for the computation of the char- 
acteristic curves, the volumes are not strictly conserved in the phase space, 
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because the scheme does not satisfy the discrete divergence free condition 



(49). 



4D Finite Volume form of the PSM scheme:. 

• PSM ID advection of in direction v\\ with velocity a^^ and 
time step At*/2 to obtain /(x, , t^H'^^). 

• PSM ID advection of f{x, f y , t^ll/^) in direction z with velocity a" and 
time step At* /2 to obtain /(x, , t^/^). 

• PSM 2D advection of /(x, , t^/^) in each plane {r,9) with velocities 
{a^,ag) and time step At* to obtain /(x, , f '^). 

• PSM ID advection of /(x, , f'^) in direction z with velocity a" and 
time step At*/2 to obtain f{x,vpt^). 

• PSM ID advection of /(x,U||,t^) in direction tiy with velocity a^^^ and 
time step At* /2 to obtain /(x,?;||,ni) = /(x, wy , 

Each PSM ID advection is achieved using the standard ID semi-Lagrangian 



scheme as described in section 2.2.2 The PSM 2D advection in (r,9) is 



achieved with the Finite Volume form as described in section 13.21 Since 



we use here the scheme 3.2, the volumes are strictly conserved in the phase 



space, because the scheme does satisfy the discrete divergence free condition 



(49). Even if we use the semi-Lagrangian PSM ID advection in z and uy 
directions, the property is kept because the velocity is constant in these 
directions. 

4-4- Results 

The mesh is 128 x 256 x 32 x 16 cells in r, 9, z, wy directions. Boundary 
conditions are periodic for directions 9 and z and Neumann [df/dn = 0) in r 
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and I'll. We first ran the reference test case witli the non-conservative Back- 



ward Semi-Lagrangian (BSL) scheme in section 2.2.1 which is currently used 
in the GYSELA code. Then we ran four test cases to show the influence of 
each numerical treatment: the standard conservative semi-Lagrangian PSM 



scheme in section 4.3.2 with ID directional splitting (PSM Directional Split- 
ting ID) and the same with the SLS limiter (SLS Directional Splitting ID), 



the unsplit Finite Volume form of the PSM scheme in section 4.3.2 (PSM 
Finite Volume) and the same with the SLS limiter (SLS Finite Volume). 
The computed 4D distribution functions are pictured in Fig. [5] at time 
t = 1800 and in Fig. [6] at time t = 4400. We only present 2D slices {r,9) of 
the distribution function at = and for a given value of z = zq. In these 
figures X stands for r direction and Y for the 6 direction. At initial time 
t = 0, the minimum and maximum values of the distribution function in 
this slice are (min. = 0.331, max. = 0.4187) and these values should be the 
same at any time of the computation if the maximum principle would be re- 
spected. In Fig. [5} we show pictures of each scheme result at time t = 1800, 
which corresponds approximately to the beginning of the non-linear tur- 
bulent phase saturation. Small structures are appearing and interact with 
each others. All results are still close qualitatively. However, we already see 
oscillations in the solution obtained with PSM DS (Directional Splitting), 
where the minimum and maximum values (min . = 0.3086, max . = 0.4427) 
are already quite different than the one at initial time. The PSM Finite Vol- 
ume (PSM FV) form and the SLS DS better keep these extrema, but only 
SLS Finite Volume keep the extrema unchanged until time t = 1800 with 
a really similar behaviour of the solution. In Fig. [6], we show pictures of 
each scheme result at time t = 4400 when turbulence is well developed. We 
see that the standard PSM scheme creates a lot of unphysical oscillations 
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Figure 5: Simulation with 128x256x32x16 cells — PSM Directional Splitting ID (up- 
left) — PSM Finite Volume (up-right) — SLS Directional Splitting ID (down-left) — 
SLS Finite Volume (down-right) — time =1800. 
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(structures are reaching the boundaries in r) and may crash the compu- 
tation. The PSM FV form and the SLS DS better keep the turbulence 
structures, but stih oscihations are created. The SLS Finite Volume keep 
the extrema of the solution reasonably well (min . = 0.3263, max . = 0.4221) 
(the SLS limiter does not provide a maximum principle) and the solution 
is smooth. We may say that the added diffusion with the limiter helps the 
scheme to diffuse subgrid structures without creating oscillations. The di- 
vergence free property of the Finite Volume scheme is important to cure to 
solution from instabilities that can be seen at r values close the average value 
of r (vertical line at the middle of pictures in Fig. [g]) in the SLS Directional 
Splitting solution compare to the SLS Finite Volume solution. In Fig. [TJ 
we see in the reference BSL solution at time t = 1800 spurious oscillations 
produced during the reconstruction step of the distribution function, which 
is the only possibility to break the maximum principle for the BSL scheme: 
here extrema are (min . = 0.3124, max . = 0.4430) instead of values at initial 
time (min. = 0.331, max. = 0.4187). At time t = 4400, we see spurious 
oscillations as well, but the maximum principle is better satisfied than with 
the standard PSM DS scheme in Fig. [6| because no conservation of volumes 



in the phase space has to be satisfied, as it is explained in section 2.2.1 



5. Conclusion and perspectives 

The PSM scheme has been successfully integrated in the GYSELA code 
and has been tested on 4D Drift-Kinetic test cases. We had first experimen- 
tally stated and afterward explained in this paper that the PSM scheme can 
be unstable without taking care of a velocity field divergence free condition. 
The numerical results show that the study of the volume evolution in the 
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Figure 6: Simulation with 128x256x32x16 cells — PSM Directional Splitting ID (up- 
left) — PSM Finite Volume (up-right) — SLS Directional Splitting ID (down-left) — 
SLS Finite Volume (down-right) — time =4400, with all color tables set to the minimum 
and maximum value at initial time. r,o 




Figure 7: Reference BSL Simulation with 128x256x32x16 cells — BSL at time=1800 
with the real color table values (left) — BSL at time=4400 with the color table set to the 
minimum and maximum value at initial time (left). 
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phase space is fruitful. Notice that this conservative scheme properly al- 
lows a directional splitting, in the semi-Lagrangian or in the Finite Volume 
form, what is not the case with the BSL scheme. The Slope Limited Splines 
(SLS) limiter is efficient to cut off spurious oscillations of the standard PSM 
scheme by adding diffusion that helps eventually the scheme to manage 
small structures below the cell size. Of course, the PSM scheme should be 
further validated as well as its integration in the GYSELA code using the 
gyrokinetic 5D model in toroidal geometry. In particular, the curvature of 
the mesh couple several directions by the geometrical Jacobian which makes 
the divergence free condition more complex, as well as the writing of the 
Quasi-Neutral solver and the Gyroaverage operator. 
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